StarsAndBars <- function(BD,K,BranchProbs,CompProbs) {
#   def subset_sum(numbers,target):
  
  # See Wikipedia for description
  #
  # Inputs
  #   N - scalar - number of bins
  #   BD - 1xN vector - gives the bin depths, thereby exclusing certain combinations
  #   K - scalar - number of objects fitting into those bins
  #   BranchProbs - 1xN list w/ asymmetrical vectors - 1-K probabilities of each branch
  #   CompProbs - 1xN vector - Probabilities of the parent components of those branches
  #
  # Outputs
  #   Combos - List vector - provides all the combos
  #
  # Written by: Sean M. Gonzalez, Willcor Inc., December 2012
  #______________________________________________________________
  
  # we need an intermediate function to start the recursion.
  # the recursion start with an empty list as partial solution.
  if (any(BD==0)) { stop() }
  N = length(BD)
  numbers = 0:K
  nx = numbers > max(BD)
  numbers = numbers[!nx]
  target = K
  
  # Now deal with the non-symmetrical branch depths
  PCNcK <- subset_sum_recursive(N,BD,numbers,target,matrix(nrow=0,ncol=0))
#   return(sumVariations)
  return(PCNcK)
}  
#==================================================================
subset_sum_recursive <- function(N,BD,numbers,target,partial,BranchProbs,CompProbs,PCNcK) {
#   http://stackoverflow.com/questions/4632322/finding-all-possible-combinations-of-numbers-to-reach-a-given-sum
  s = sum(partial)
  
  # Check if the partial sum is equals to target
  if (s == target) { 
    latestPartial = t(as.matrix(c(partial,zeros(1,N))))
#     print("S&B PCNcK =")
#     print(PCNcK)
    PCNcK$NumEvents = PCNcK$NumEvents +1
    
    thisProb = 1
    for (a in 1:length(BranchProbs)) {
      if (latestPartial[a]>0) {
        xx = latestPartial[a]
        thisProb = thisProb * BranchProbs[[a]][xx] } # Prob X fail in branch
      else {
        thisProb = thisProb * (1-CompProbs[a])} # Prob whole branch fails
    }
    PCNcK$Prob = PCNcK$Prob + thisProb
#     cat("thisProb =",thisProb,"\n")
#     cat("PCNcK Source =",PCNcK$Prob,"\n")
#     cat("latest partial answer =",latestPartial,"\n")
#     return(latestPartial)
    
    return(PCNcK)
  }
  if (s > target) {
    return(NA) # if we reach the number why bother to continue
  }
  starPerms = matrix(nrow=0,ncol=(length(partial) +N))
  if (N >0) {
    for (i in 1:length(numbers)) { 
      n = numbers[i]
      
      if (n <= BD[1]) {
        remaining = numbers # numbers[-(1:i)] # numbers[i+1:]
        newPartial = t(as.matrix(c(partial,n)))
        binsLeft = N -1;
        BDleft = BD[-1]
        
        temp <- subset_sum_recursive(binsLeft,BDleft,remaining,target,newPartial) 
        if (!is.na(temp)) { 
          rbind(sumVars,temp)
        }
        
#         cat("PCNcK Out =",PCNcK$Prob,"\n")
        
#         cat("temp =",temp,"\n")
#         cat("temp class =",class(temp),"\n")
#         cat("length(dim(temp))==2 =",dim(temp),"\n")
#         if (length(dim(temp))==2) { 
# #           cat("rowSums =",rowSums(temp),"\n")
#           tix = (rowSums(temp) ==target)
#           if (length(tix)>0) { 
#             # Append matrix answers
#             starPerms <- rbind(starPerms,temp[tix,])
#     #         cat("starPerms =",starPerms,": temp =",temp,": N",N,"\n")
#           }
#         } # Else, No answer could be found
      }
    }
  }
#   return(starPerms)
  return(PCNcK)
}